Using time series vector features for annual cultivated land mapping: A trial in northern Henan, China

Annual monitoring of the spatial distribution of cultivated land is important for maintaining the ecological environment, achieving a status quo of land resource management, and guaranteeing agricultural production. With the gradual development of remote sensing technology, it has become a common practice to obtain cultivated land boundary information on a large scale with the help of satellite Earth observation images. Traditional land use classification methods are affected by multiple types of land cover, which leads to a decrease in the accuracy of cultivated land mapping. In contrast, although the current advanced methods (such as deep learning) can obtain more accurate cultivated land mapping results than traditional methods, such methods often require the use of a massive amount of training samples, large computing power, and highly complex model tuning processes, increasing the cost of mapping and requiring the involvement of more professionals. This has hindered the promotion of related methods in mapping institutions. This paper proposes a method based on time series vector features (MTVF), which uses vector thinking to establish the features. The advantage of this method is that the introduction of vector features enlarges the differences between the different land cover types, which overcomes the loss of mapping accuracy caused by the influences of the spectra of different ground objects and ensures the calculation efficiency. Moreover, the MTVF uses a traditional method (random forest) as the classification core, which makes the MTVF less demanding than advanced methods in terms of the number of training samples. Sentinel-2 satellite images were used to carry out cultivated land mapping for 2020 in northern Henan Province, China. The results show that the MTVF has the potential to accurately identify cultivated land. Furthermore, the overall accuracy, producer accuracy, and user accuracy of the overall study area and four sub-study areas were all greater than 90%. In addition, the cultivated land mapping accuracy of the MTVF is significantly better than that of the maximum likelihood, support vector machine, and artificial neural network methods.

Introduction Cultivated land is one of the main types of land cover, and it is a key component of human food production. At present, the world's cultivated land feeds more than seven billion people. Due to rapid population growth, the development of cultivated land far exceeds its carrying capacity. Therefore, monitoring the annual spatial distribution of cultivated land is an important prerequisite for protecting the ecological environment [1], establishing a new status quo in land resource management [2], and guaranteeing agricultural production [3][4][5][6][7][8]. The mastery of annual cultivated land mapping data has become an important research topic.
Satellite remote sensing technology has gradually become one of the main methods for cultivated land mapping due to its timeliness, low cost, and large-scale observation capabilities. At present, with the help of remote sensing observation archives, many global cultivated land mapping products and global land cover mapping products containing cultivated land map layers have been released, including the finer resolution observation and monitoring of global land cover (FROM-GLC) [9,10], GlobeLand 30 [11], data and information system global land cover (DISCover) [12], and moderate resolution imaging spectroradiometer (MODIS) [13,14] land cover products. However, the above products are only available for certain years (such as 2015 and 2020), and thus, they do not meet the requirements for annual cultivated land mapping. Furthermore, due to the limitations of the spatial scales of the existing global cultivated land mapping products, the definition of cultivated land often fails to take into account the local scale. Sustainable development research, food security, and other fields have created new requirements for higher resolution, high precision, and local cultivated land monitoring data.
There are various methods for cultivated land mapping based on remote sensing observation data, including traditional supervised classification models (e.g., the maximum likelihood method [15], support vector machines [16], and artificial neural networks [17]) and advanced methods (e.g., deep learning [18] and artificial immune networks [19]). However, regarding the cultivated land mapping process, the existing methods still encounter many limitations. For example, traditional methods are susceptible to the effects of multiple land cover types, resulting in a decrease in the accuracy of cultivated land surveying and mapping [20]; whereas advanced methods (e.g., deep learning) need to rely more on massive training samples, supercomputing power, and complex model tuning [21,22]. Establishing a simple but reliable method is the key to cultivated land mapping research.
The selection of data for cultivated land mapping is also important. One feasible idea is to use time series data. Many studies have shown that multi-spectral remote sensing images based on time series sequences are an effective means of large-scale, long-term, continuous agricultural remote sensing mapping [23,24]. Multispectral remote sensing image data based on time series can overcome the influences of various factors, such as the weather, and can provide a data basis for continuous crop growth curve extraction [25,26]. Time series analysis combined with the vegetation index is also an effective idea for cultivated land mapping. As the most widely used characteristic parameter to describe vegetation phenological changes [27], the vegetation index time series can reflect the dynamic changes in different crop types over time. The vegetation index based on time series multi-spectral remote sensing image data reflects the dynamic changes in different crop types over time.
Based on the above discussion, a method based on time series vector features (MTVF) is proposed in this article. The core of this method is to use a time series based on spectral and vegetation indices as a vector and to extract the vector features to distinguish the differences between cultivated land and other land cover types. The purposes of this approach are as follows. 1) This approach can widen the difference between other land cover types and cultivated land and reduce the impact of spectral overlap. 2) Compared with deep learning methods, the proposed method requires fewer training samples and less computing power, and thus, it can significantly improve the efficiency of cultivated land mapping. In this study, Sentinel-2 satellite data were used to generate a 10-m spatial resolution cultivated land map because these data have a short revisit period and contain rich red-edge band information. This method was applied to the study area in the northern part of Henan Province, China, to verify its applicability to cultivated land mapping. In addition, four sub-study areas located in different landscapes within the study area were established to evaluate the cultivated land maps. Finally, the MTVF was compared with three traditional supervised classification models (the maximum likelihood, support vector machine, and artificial neural network models) to evaluate the advantages of this method in cultivated land mapping.

Study area
In this study, a 50 km × 50 km area in the northern part of Henan Province was selected as the study area (Fig 1). The selection of this location was mainly based on the following factors: (i) the diversity of the plant types, (ii) the complexity of the land cover, (iii) the presence of multiple types of topography, and (iv) the heterogeneous spatial distribution of the grain yield.
Henan Province, which is located in the central part of China, is one of the country's major grain cultivation areas. According to official statistics [28], the total cultivated land area in Henan Province was about 14.68×10 6 hectares in 2019, an increase of about 2.01% compared to that in 2012. The northern part of Henan Province has a warm, temperate, monsoon climate, with an average annual rainfall of 500-900 mm and an average annual sunshine duration of 1285.7-2292.9 h, which makes it suitable for the growth of a variety of crops. Agricultural production activities in this area are affected by many factors, such as the water resources [29], labor per unit area, and urbanization rate [30].
To evaluate the effectiveness of the cultivated land mapping, four sub-areas within the study area were selected. Sub-study area A (LZ) was located in Linzhou County. LZ was located in the Taihang Mountains, and it contained a large amount of abandoned land, as well as abundant grass, woodland, and other natural vegetation. This was the main reason that LZ was selected as a sub-study area. Sub-study area B (XX) was located in the city of Xinxiang. The reason for choosing XX was mainly to consider the impact of open-pit mines on cultivated land. Sub-study area C (HB) was located in the city of Hebi. The natural grassland in HB was very lush, and the distribution of the cultivated land landscape was relatively fragmented. Experimentation in this area helped us to analyze the ability of the algorithm to separate grassland and cultivated land. Sub-study area D (WH) was located in Weihui County and contained large grape plantations. WH was selected to analyze the ability of the algorithm to strip orchards (grapes) in the process of cultivated land mapping. All of the sub-study areas were delimited by a 3 km × 3 km rectangle, and the winter wheat-summer corn rotation pattern was dominant in the four sub-study areas.

Satellite data
Sentinel-2 is a polar-orbit, high-resolution, multi-spectral imaging mission for terrestrial monitoring. It consists of two satellites, Sentinel-2A and Sentinel-2B, which are each equipped with a multi-spectral imager (MSI). Sentinel-2 satellite data includes data for 13 spectral bands, ranging from visible and near-infrared light to short-wave infrared light, with a ground resolution of up to 10 m (Table 1). Since the Sentinel-2 satellite can provide data in three spectral bands within the red-edge range and the data update cycle is 5 days (two satellites for monitoring), the Sentinel-2 satellite data have advantages relevant to vegetation time series monitoring.
All of the data collected by the Sentinel-2 satellite can be downloaded from the European Space Agency (ESA) Copernicus Open Access Hub (https://scihub.copernicus.eu/dhus/ #/home).
Since the main planting pattern in the study area was the rotation of winter wheat and summer corn, October (the winter wheat was sown in the study area during this month) was selected as the annual time node. Based on this idea, we selected 12 scenes of Sentinel-2 images from October 2019 to September 2020 ( Table 2). The scenes were selected to be as cloudless, fogless, and evenly distributed in each month as possible. These image data served as the basis for constructing the vegetation index time series.
It should be noted that Sentinel-2 images acquired in the middle of each month were selected to ensure the uniform distribution of the time series. Still, the acquisition times of some of the images changed to avoid excessive cloud cover. For example, the acquisition times of the October and December 2019 images were later in the month, and the acquisition times of the June and July 2020 images were in the early part of the month. The L2A level Sentinel-2 images were obtained from the ESA Copernicus Open Access Hub and had already been subjected to official ESA geometric and atmospheric corrections (see https://sentinel.esa.int/web/ sentinel/user-guides/sentinel-2-msi/processing-levels/level -2), so they could be directly used to calculate the vegetation index.

Collection of reference data
Reference land cover classes in study area. To investigate the algorithm's ability to separate cultivated land from other land cover types, based on the results of field investigations and the land cover classification system developed by the Chinese Academy of Sciences, the main land cover types in the study area were determined.
• Cultivated land: This refers to the land on which crops are grown, including food crops and some vegetables. It should be noted that fruit plantations were not included in this category, because in the study area, orchards were not the main form of agriculture and were also under the jurisdiction of the local forestry management department. • Woodland: This refers to forestry land where trees, shrubs, bamboo, and orchards grow.
• Grassland: This refers to land dominated by natural herbaceous plants. The study area did not contain artificial grassland, so it was not considered.
• Water: This refers to natural terrestrial water area and land used for water conservancy facilities, such as rivers, artificial canals, lakes, ponds, and pools.
• Artificial construction land: This refers to urban, rural, industrial, mining, and residential land.
• Bare land: This refers to land and rocks with vegetation coverage of less than 5%. In the study area, rock exposed by mining was also included in this type.

Collection of training and validation samples.
We collected training and validation samples in the study area in the form of pixels. In order to obtain the above information, we combined a field survey (May 2021) and a visual inspection of very-high-resolution (VHR) images from Google Earth and unmanned aerial vehicle (UAV) field measurements, the VHR image from Google Earth was acquired on December 31, 2019, and the VHR images from the UAV were acquired on May 11-14, 2021. December and May were chosen because of the agricultural characteristics of the study area. December is the emergence period of winter wheat, so it is possible to effectively distinguish cultivated land from other vegetation in this month. May is the maturity period of wheat, so it is convenient to screen non-vegetation covered land types in this month. During the field investigations, the coordinates and attribute information about the ground objects were mainly recorded using a handheld global positioning system (GPS). To overcome the inaccuracy caused by the different years, the samples were screened by manually comparing the 2021 VHR images acquired using the UAV and the Sentinel-2 images to ensure that the attributes of the sample points were consistent with the actual land cover types in 2020. In the study area, 698 and 14,158 well disturbed pixels were used as training samples and verification samples, respectively (Fig 2). Although including more training samples would have had a positive effect on the remote sensing supervised classification calculations, we used fewer training samples in order to verify the advantages of our algorithm. For supervised classification, the general rule is that the number of training samples per class should be 10-30 times the input image band [31][32][33]. Therefore, in this paper, ten times the number of bands of the Sentinel-2 satellite imagery was used as the reference value. Training samples for each surface coverage category were selected from the collected samples, and the remaining samples were used as the verification samples. However, the distributions of the water and bare land in the study area were limited, so only 70 pixels were selected as training samples for each of these two land cover types.
In order to verify the accuracy of the results for the four sub-study areas, several pixels were collected from each sub-study area as verification samples. According to the ground truth, these pixels were labeled as 1 (cultivated land) or 0 (non-cultivated land). The distribution of the samples is shown in Table 3.  (Table 4) were used as vectors for the subsequent calculations (see Vector Construction). Five parameters (i.e., the cosine, distance, maximum, minimum, and range of each vector) were calculated (see Feature Extraction). Then, we calculated the importance score of each parameter (based on the random forest model) to evaluate the contribution of each parameter to the cultivated land mapping and to provide a thorough analysis and discussion. By sorting and grouping the importance scores of the parameters and calculating the performances of the different groupings in terms of the random forest  prediction accuracy and the out-of-bag error, several optimal parameters that could be used in the calculations for the cultivated land mapping in the study area were determined. The number of trees was also determined in this manner (see Random Forest Model). The training samples and optimal parameters were introduced into the random forest classification algorithm to generate the cultivated land map layer in the next step. Finally, the classification accuracy of the cultivated land map was evaluated through the verification sample and compared with the results obtained using other algorithms (i.e., the maximum likelihood, support vector machine, and artificial neural network methods; see Other Classification Models and Accuracy Assessment).

Vector construction
The parameters used to build the vector included the Sentinel-2 satellite band value and vegetation indices. The reflectivity characteristics of the different features in each wavelength range were dissimilar. Thus, the 13 bands of the Sentinel-2 satellite were used as the radiation parameters to establish the vector, and these data were normalized. Furthermore, a total of 23 vegetation indices (Table 4) were selected as the vegetation index parameters for the cultivated land mapping. The selection of the vegetation indices mainly followed two rules: 1.) they can be produced from the bands of the Sentinel-2 satellite images; and 2.) they have been verified in practice. These vegetation indices were calculated using the ESA Sentinel Application Platform (SNAP) software (http://step.esa.int/main/toolboxes/snap/). All of the radiation parameters and vegetation index parameters were constructed as vectors in the form of time series and can be expressed as follows: where V p ! is the time series vector of parameter p from October 2019 to September 2020.

Feature extraction
Vector features can be used as significant parameters for remote sensing image classification, including the vector angle [56,57], vector distance [58], and extreme value [59][60][61]. Therefore, we established five parameters: the cosine, distance, maximum, minimum, and range based on the vector V p ! .

Random forest model
The random forest (RF) model is a machine learning classifier that combines multiple decision trees [62]. The random forest classification process can be described as follows: Step 1. The bootstrap sampling method is used to randomly select training samples with replacement.
Step 2. Set the corresponding decision tree model for each training sample and continue to split until all the training samples of the node are of the same type.
Step 3. The generated multiple decision trees are formed into a random forest, and the optimal classification is determined according to the voting probability.
In this study, the importance score, out-of-bag error, prediction accuracy, and classification calculations were performed using the random forest model. The importance score was used to determine the contribution of each feature to the cultivated land mapping. The out-of-bag error and prediction accuracy were used to select the feature set and to determine the number of trees in the random forest. All of the calculations were based on the Scikit-learn [63] machine learning algorithm library built in the Python language. Importance score. The Gini coefficient generated by the random forest algorithm was used to compare the contributions of the individual features to the cultivated land mapping. The classification and regression tree method was integrated with the RF to provide the Gini coefficient for the next split. Thus, the importance of each feature was expressed as follows: where G is the Gini coefficient before the split, and G L and G R are the Gini coefficients of the left and right branches after the split. It is assumed that there are a total of n trees in the forest that use feature f, and the number of splits of each tree is k. I f is the importance, and I is the sum of the importance of all of the features. S f is the importance score of feature f, which is the normalized value of I f . Out-of-bag error. A significant advantage of the random forest model is that it can build an unbiased estimate of the error internally, which is called the out-of-bag error (oob). For each tree, the randomly selected samples were approximately 63.2% of the total number of samples, and the remaining approximately 36.8% of the samples were designated as the oob samples of this tree. The oob of a random forest is the mean value of the oob of all of the trees in the model.

Feature selection
An excessive number of features leads to higher computational costs and redundancy, so it is necessary to select appropriate features to participate in the next step of the calculation. In this study, the feature selection mainly included three steps. First, all of the parameters were sorted according to their importance scores to obtain their distribution characteristics. Then, all of the parameters were grouped according to their distribution characteristics to obtain several feature groups. Finally, all of the feature groups were input into the random forest model, the performances of the prediction accuracy and oob were calculated for different numbers of trees, and the group with the highest prediction accuracy and the smallest out-of-band error was selected as the feature group.

Other classification models
To evaluate the cultivated land mapping model developed in this study, the maximum likelihood model, support vector machine model, and artificial neural network model were introduced for comparison. The same training and validation samples used in the MTVF were used in all of the traditional models to ensure the objectivity of the cultivated land mapping results of the comparison of the different models. The values of the parameters of each model were set within the range recommended by the developer because this is generally considered to be a safe practice [64].
Maximum likelihood model. The maximum likelihood model is one of the most widely used supervised classification models. It is a theoretical point estimation algorithm. The maximum likelihood model assumes that the distributions of the various types of data in each band are Gaussian, and each peak represents a unique feature category. Based on the training samples, the statistics of the normal distribution are obtained, and then, the probability of each pixel belonging to a different normal distribution is calculated. Finally, the pixel is classified into the category with the largest probability. The classification results of the maximum likelihood model have the advantages of stability, reliability, algorithm simplicity, high accuracy, and fast calculation times [15]. Support vector machine model. The concept of the support vector machine (SVM) was first proposed by Cortes and Vapnik [16]. It maps the vector to a higher dimension and implements classification by introducing a kernel function to map the sample data in the lowdimensional feature space to the high-dimensional feature space. In recent years, support vector machine models have been widely used for the segmentation, fusion, and classification of high-spatial-resolution images [65].
Artificial neural network model. The artificial neural network model is a non-parametric classification method with a good adaptability and a complex mapping capability. By mimicking the brain's structure and functions, it realizes non-linear data pattern recognition and can effectively combine the spectral and textural features of images to improve the classification accuracy [17,66,67].

Accuracy assessment
A confusion matrix was used to assess the accuracy of the cultivated land mapping. In this study, the accuracy of the producer (PA), the accuracy of the user (UA), the overall accuracy (OA), the kappa coefficient of variation (kappa), and the ground truth (by pixel) were chosen as the indices to measure the cultivated land mapping ability of each algorithm.

Importance score of features
By calculating a total of 180 features, the 50 most important features are listed in Table 5. Among all of the features, Band 12 and Band 8A exhibited advantages in cultivated land mapping, with importance scores of 0.1723 and 0.1368, respectively. Among the top 10 parameters (importance scores of > 0.02), five of the parameters were based on the band value of the image, and the other five parameters were based on the vegetation indexes. The parameters characterized by the vector angle and the vector distance occupied dominant positions in the top 10 ranking of the importance. For the vegetation indexes, the normalized difference vegetation index (NDVI 705 , importance score of 0.0526) made the largest contribution to the cultivated land mapping, followed by the atmospherically resistant vegetation index (ARVI, importance score of 0.0516) and green normalized difference vegetation index (GNDVI, importance score of 0.0440). There were 26 features with importance scores greater than the average score (1/180 � 0.0056).

Feature selection results
The visualization of the importance score from high to low (Fig 4) shows that the importance scores exhibited three nodes: 0.1, 0.04, and 0.02. Based on this, four feature groups were established in this study: Group A (importance scores of > 0.1), Group B (importance scores of > 0.04), Group C (importance scores of > 0.02), and Group D (importance scores of > 0). The four feature groups were input into the random forest model to obtain the prediction accuracy and to determine the changes in the oob (Fig 5). The results revealed that Group C had the highest prediction accuracy and the lowest oob. Therefore, the features with importance scores of > 0.02 were selected, and the number of trees was set to 400.

Accuracy assessment of cultivated land mapping
Entire study area. The MTVF exhibited a strong ability to map cultivated land in the study area, and its OA, PA, UA, and kappa coefficient performance reached higher levels (accuracy > 90%, kappa coefficient > 0.85) ( Table 6). Compared with other land types, the mapping accuracies (PA and UA) of the cultivated land were both the highest. The classification accuracies (PA and UA) of the grassland, woodland, and water bodies were also relatively high. The classification accuracies (PA and UA) of the artificial construction land and bare land were worse than those of the other land cover types. In particular, the UA of the artificial construction land was 73.95%, and the PA of the bare land was 61.67%. Thus, the classification of these land cover types was considered to have a poor accuracy (accuracy <75%). These land cover types may have been affected by the input features because the input features were selected specifically for cultivated land mapping.
Sub-study areas. In the four sub-study areas, the MTVF also yielded a strong cultivated land mapping accuracy ( Table 7). The mapping accuracies of the cultivated land in the LZ and HB areas were the highest of the four sub-study areas, indicating the robustness of the method in distinguishing cultivated land from grassland and woodland. The UA and PA values of area XX were the lowest of the four sub-study areas, indicating that some of the cultivated land was incorrectly classified into other ground feature categories. The reason for the lower accuracy in area XX was that this area was also affected by dense and fragmented patterns and mixed

PLOS ONE
pixels. The results for area WH show that the method proposed in this paper is affected by artificial orchard features to a certain extent.
Comparison with other classification models Fig 6 shows the difference between the cultivated land mapping results of the MTVF and the other models (i.e., the maximum likelihood (ML), support vector machine, and artificial neural network (ANN) models). From the perspective of the entire study area, the MTVF and ML had the best cultivated land mapping effects, and the large area of cultivated land located in a flat area could be distinguished with distinct boundaries. The difference between the cultivated land mapping capabilities of the MTVF and ML was mainly reflected in the mountainous areas (northwestern part of the entire study area), in which the ML contained more speckle noise, and much of the grassland was misclassified as cultivated land. In contrast, the ANN and SVM could not effectively filter out the cultivated land pixels in the study area, and the SVM had the worst cultivated land mapping ability. In the four sub-study areas, the SVM had the worst cultivated land mapping ability of all of the models and was not included in the subsequent comparison. In area LZ, the MTVF effectively distinguished several regular artificial forest land areas connected to cultivated land. Although the ANN and ML also exhibited this ability, some of the plantation pixels were still mistakenly classified as cultivated land. In area XX, broken patches and mixed pixels were the main factors that affected the accuracy of the cultivated land mapping. The ANN and ML had a limited ability to overcome these influencing factors. Therefore, the ANN and ML misclassified woodland and grassland in this substudy area and generated more speckle noise. In area HB, the MTVF, ML, and ANN effectively avoided the influence of the grassland, but only the MTVF distinguished between the woodland pixels and cultivated land pixels. In the WH, the MTVF effectively avoided the impact of the roads between the cultivated land and retained wider country roads. However, the MTVF failed to completely avoid the impact of the orchards and the orchards were mistakenly classified as cultivated land. The MTVF was compared with other classification models. From the perspective of the classification accuracy (Table 8), the MTVF, ANN, and ML all produced cultivated land maps with good accuracies. Their accuracies were greater than 75%, and their kappa coefficients were greater than 0.75. The cultivated land map generated by the SVM model had the worst accuracy. The SVM's accuracy was not only the worst for the entire study area, but its accuracy performances in the four sub-study areas were also inferior to those of the other methods. In sub-study area XX, although the ANN achieved the highest OA and kappa coefficients, the MTVF achieved better UA and PA values. In general, the MTVF exhibited a better performance in cultivated land mapping than the other models.
In this study, the statistics of the ground truth (pixels) in the confusion matrices of the cultivated land maps produced using the various methods were computed ( Table 9). The cultivated land maps produced using the MTVF were the least affected by the other types of ground features. The impact of the woodland was dominant, but the artificial construction land and grassland also had impacts. Factors such as artificially planted forests in cultivated land, field ravines, and rural residential land boundaries were the reasons for the loss of accuracy of the MTVF. The ANN and SVM exhibited a weak ability to avoid the influence of the woodland, and the ANN model was also more susceptible to the influence of the grassland. In general, the forest land, grassland, and artificial construction land were the factors that affected the mapping of the cultivated land because there were often a large number of mixed pixels between the cultivated land in northern Henan and the above surface cover types, which led to misclassification. By limiting the impacts of these types of land, the MTVF achieved better results.

Discussion
In this study, vector thinking was a core concept, and an MTVF that uses time series multispectral data for annual cultivated land mapping was developed. The feature extraction of the time series vectors was the core of this research. In this study, by introducing unit vectors, we extracted five features of each time series vector: Max p (the maximum value inside vector p), Min p (the minimum value inside vector p), Ran p (the difference between Max p and Min p ), Cos p (the cosine of the angle between time series vector p and the unit vector), and Dis p (the distance between time series vector p and the unit vector). The features extracted from the different time series vectors effectively reduced the misclassification of land cover types caused by similar spectra and enlarged the differences between the land cover types and cultivated land. In the MTVF, the selection of the features is very important. Selecting too few features cannot provide a sufficient classification basis for the classification model, and selecting too many features will lead to redundancy. In this study, the Gini coefficient was introduced to measure the importance of the different characteristics to the cultivated land mapping. Statistical analysis of the importance of each feature revealed that B12_Cos (importance score of 0.1723) and B8A_Dis (importance score of 0.1368) had the most significant Gini importance scores, indicating that the shortwave infrared band (Band 12) and the red edge band (Band 8A) made relatively large contributions to the cultivated land mapping. Among the vegetation index parameters, the NDVI 705 _Min (importance score of 0.0526, ranked third) and ARVI_Cos (importance score of 0.0516, ranked fourth) had the highest importance scores. Traditional vegetation indices such as the normalized vegetation index (NDVI) and the enhanced vegetation index (EVI) did not have high importance scores in this study. This phenomenon revealed that non-traditional vegetation indices can also be valuable in cultivated land mapping, which is similar to the results of several previous studies [68,69]. The parameters characterized by the vector angle and vector distance dominated the 10 most important parameters, indicating that spatial features based on time series vectors can reflect the differences between cultivated land and other land use types better than vector extreme value features. It should be noted that for different research areas, changes in the types of land cover will also lead to changes in the importance of the features. When conducting cultivated land mapping in different research areas, it is necessary to reassess the importance of these features [70].
In this study, the accuracy of the cultivated land mapping conducted using the MTVF was evaluated in the entire study area and in four sub-study areas. It was found that the cultivated land map obtained using the MTVF had the highest accuracy compared to the traditional classification models (i.e., the maximum likelihood, artificial neural network, and support vector machine models). The cultivated land map obtained using the MTVF had the lowest error and limited the effects of mixed pixels to a certain extent (Fig 6). However, the MTVF was still subject to interference from mixed pixels, resulting in some uncertainty in its ability to perform cultivate land mapping. In the study area, affected by the traditional farming practices in northern China, the distribution of cultivated land in some areas was irregular. There was cultivated land within villages and near mountain ravines, and trees were also planted in the cultivated land. The aforementioned areas with heterogeneous landscapes led to a large number of mixed pixels, resulting in errors in the mapping of the cultivated land [19,71]. The spectral confusion between the woodland and cultivated land also contributed to the uncertainty of the cultivated land map obtained using the MTVF. In sub-study area WH, the spectral confusion between the cultivated land and orchards contributed to the misclassification of the MTVF to some extent (Fig 6).

Conclusions
In this study, a simple and effective method for cultivated land mapping was developed. The MTVF has a stronger ability to eliminate the influences of other vegetation. By introducing vector thinking, an MTVF was developed. The cultivated land mapping performance of the MTVF was evaluated in the entire study area and in four sub-study areas, The main conclusions of this study are as follows.
1. The MTVF has a high potential for cultivated land mapping and achieved a high accuracy (greater than 90%) in the study area. The MTVF mainly mixed the effects of the pixels in the mapping of the cultivated land, especially where the land cover was complicated. However, in some cases, the MTVF also had the ability to limit the influence of the mixed pixels.
2. In terms of the importance scores, the shortwave infrared and red-edge bands of the Sentinel-2 satellite have a high potential for cultivated land mapping. The non-traditional vegetation indices were superior to the traditional vegetation indices in terms of cultivated land mapping. The spatial features based on time series vectors reflected the differences between cultivated land and other land use types better than the vector extreme features.
3. Compared with other models (i.e., the maximum likelihood, support vector machine, and artificial neural network models), the MTVF achieved the best results in the study area, but it still suffers interference from artificial woodland, field gullies, and rural settlement boundaries, which decrease the accuracy of the cultivated land map.
The MTVF provides a new method for cultivated land mapping. Future research should focus on combining this method with the mixed pixel decomposition algorithm. Combining multi-source sensor data (e.g., synthetic aperture radar) for use in cultivated land mapping should also be a future research focus because the mapping accuracy of the MTVF can still be affected by other ground cover types.